{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# State space modeling: Local Linear Trends"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook describes how to extend the statsmodels statespace classes to create and estimate a custom model. Here we develop a local linear trend model.\n",
    "\n",
    "The Local Linear Trend model has the form (see Durbin and Koopman 2012, Chapter 3.2 for all notation and details):\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "y_t & = \\mu_t + \\varepsilon_t \\qquad & \\varepsilon_t \\sim\n",
    "    N(0, \\sigma_\\varepsilon^2) \\\\\n",
    "\\mu_{t+1} & = \\mu_t + \\nu_t + \\xi_t & \\xi_t \\sim N(0, \\sigma_\\xi^2) \\\\\n",
    "\\nu_{t+1} & = \\nu_t + \\zeta_t & \\zeta_t \\sim N(0, \\sigma_\\zeta^2)\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "It is easy to see that this can be cast into state space form as:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "y_t & = \\begin{pmatrix} 1 & 0 \\end{pmatrix} \\begin{pmatrix} \\mu_t \\\\ \\nu_t \\end{pmatrix} + \\varepsilon_t \\\\\n",
    "\\begin{pmatrix} \\mu_{t+1} \\\\ \\nu_{t+1} \\end{pmatrix} & = \\begin{bmatrix} 1 & 1 \\\\ 0 & 1 \\end{bmatrix} \\begin{pmatrix} \\mu_t \\\\ \\nu_t \\end{pmatrix} + \\begin{pmatrix} \\xi_t \\\\ \\zeta_t \\end{pmatrix}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Notice that much of the state space representation is composed of known values; in fact the only parts in which parameters to be estimated appear are in the variance / covariance matrices:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "H_t & = \\begin{bmatrix} \\sigma_\\varepsilon^2 \\end{bmatrix} \\\\\n",
    "Q_t & = \\begin{bmatrix} \\sigma_\\xi^2 & 0 \\\\ 0 & \\sigma_\\zeta^2 \\end{bmatrix}\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.stats import norm\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To take advantage of the existing infrastructure, including Kalman filtering and maximum likelihood estimation, we create a new class which extends from `statsmodels.tsa.statespace.MLEModel`. There are a number of things that must be specified:\n",
    "\n",
    "1. **k_states**, **k_posdef**: These two parameters must be provided to the base classes in initialization. The inform the statespace model about the size of, respectively, the state vector, above $\\begin{pmatrix} \\mu_t & \\nu_t \\end{pmatrix}'$, and   the state error vector, above $\\begin{pmatrix} \\xi_t & \\zeta_t \\end{pmatrix}'$. Note that the dimension of the endogenous vector does not have to be specified, since it can be inferred from the `endog` array.\n",
    "2. **update**: The method `update`, with argument `params`, must be specified (it is used when `fit()` is called to calculate the MLE). It takes the parameters and fills them into the appropriate state space matrices. For example, below, the `params` vector contains variance parameters $\\begin{pmatrix} \\sigma_\\varepsilon^2 & \\sigma_\\xi^2 & \\sigma_\\zeta^2\\end{pmatrix}$, and the `update` method must place them in the observation and state covariance matrices. More generally, the parameter vector might be mapped into many different places in all of the statespace matrices.\n",
    "3. **statespace matrices**: by default, all state space matrices (`obs_intercept, design, obs_cov, state_intercept, transition, selection, state_cov`) are set to zeros. Values that are fixed (like the ones in the design and transition matrices here) can be set in initialization, whereas values that vary with the parameters should be set in the `update` method. Note that it is easy to forget to set the selection matrix, which is often just the identity matrix (as it is here), but not setting it will lead to a very different model (one where there is not a stochastic component to the transition equation).\n",
    "4. **start params**: start parameters must be set, even if it is just a vector of zeros, although often good start parameters can be found from the data. Maximum likelihood estimation by gradient methods (as employed here) can be sensitive to the starting parameters, so it is important to select good ones if possible. Here it does not matter too much (although as variances, they should't be set zero).\n",
    "5. **initialization**: in addition to defined state space matrices, all state space models must be initialized with the mean and variance for the initial distribution of the state vector. If the distribution is known, `initialize_known(initial_state, initial_state_cov)` can be called, or if the model is stationary (e.g. an ARMA model), `initialize_stationary` can be used. Otherwise, `initialize_approximate_diffuse` is a reasonable generic initialization (exact diffuse initialization is not yet available). Since the local linear trend model is not stationary (it is composed of random walks) and since the distribution is not generally known, we use `initialize_approximate_diffuse` below.\n",
    "\n",
    "The above are the minimum necessary for a successful model. There are also a number of things that do not have to be set, but which may be helpful or important for some applications:\n",
    "\n",
    "1. **transform / untransform**: when `fit` is called, the optimizer in the background will use gradient methods to select the parameters that maximize the likelihood function. By default it uses unbounded optimization, which means that it may select any parameter value. In many cases, that is not the desired behavior; variances, for example, cannot be negative. To get around this, the `transform` method takes the unconstrained vector of parameters provided by the optimizer and returns a constrained vector of parameters used in likelihood evaluation. `untransform` provides the reverse operation.\n",
    "2. **param_names**: this internal method can be used to set names for the estimated parameters so that e.g. the summary provides meaningful names. If not present, parameters are named `param0`, `param1`, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Univariate Local Linear Trend Model\n",
    "\"\"\"\n",
    "class LocalLinearTrend(sm.tsa.statespace.MLEModel):\n",
    "    def __init__(self, endog):\n",
    "        # Model order\n",
    "        k_states = k_posdef = 2\n",
    "\n",
    "        # Initialize the statespace\n",
    "        super(LocalLinearTrend, self).__init__(\n",
    "            endog, k_states=k_states, k_posdef=k_posdef,\n",
    "            initialization='approximate_diffuse',\n",
    "            loglikelihood_burn=k_states\n",
    "        )\n",
    "\n",
    "        # Initialize the matrices\n",
    "        self.ssm['design'] = np.array([1, 0])\n",
    "        self.ssm['transition'] = np.array([[1, 1],\n",
    "                                       [0, 1]])\n",
    "        self.ssm['selection'] = np.eye(k_states)\n",
    "\n",
    "        # Cache some indices\n",
    "        self._state_cov_idx = ('state_cov',) + np.diag_indices(k_posdef)\n",
    "\n",
    "    @property\n",
    "    def param_names(self):\n",
    "        return ['sigma2.measurement', 'sigma2.level', 'sigma2.trend']\n",
    "\n",
    "    @property\n",
    "    def start_params(self):\n",
    "        return [np.std(self.endog)]*3\n",
    "\n",
    "    def transform_params(self, unconstrained):\n",
    "        return unconstrained**2\n",
    "\n",
    "    def untransform_params(self, constrained):\n",
    "        return constrained**0.5\n",
    "\n",
    "    def update(self, params, *args, **kwargs):\n",
    "        params = super(LocalLinearTrend, self).update(params, *args, **kwargs)\n",
    "        \n",
    "        # Observation covariance\n",
    "        self.ssm['obs_cov',0,0] = params[0]\n",
    "\n",
    "        # State covariance\n",
    "        self.ssm[self._state_cov_idx] = params[1:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this simple model, we can estimate the parameters from a local linear trend model. The following example is from Commandeur and Koopman (2007), section 3.4., modeling motor vehicle fatalities in Finland."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import requests\n",
    "from io import BytesIO\n",
    "from zipfile import ZipFile\n",
    "    \n",
    "# Download the dataset\n",
    "ck = requests.get('http://staff.feweb.vu.nl/koopman/projects/ckbook/OxCodeAll.zip').content\n",
    "zipped = ZipFile(BytesIO(ck))\n",
    "df = pd.read_table(\n",
    "    BytesIO(zipped.read('OxCodeIntroStateSpaceBook/Chapter_2/NorwayFinland.txt')),\n",
    "    skiprows=1, header=None, sep='\\s+', engine='python',\n",
    "    names=['date','nf', 'ff']\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we defined the local linear trend model as extending from `MLEModel`, the `fit()` method is immediately available, just as in other statsmodels maximum likelihood classes. Similarly, the returned results class supports many of the same post-estimation results, like the `summary` method.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                           Statespace Model Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                    lff   No. Observations:                   34\n",
      "Model:               LocalLinearTrend   Log Likelihood                  27.510\n",
      "Date:                Tue, 24 Dec 2019   AIC                            -49.020\n",
      "Time:                        15:04:44   BIC                            -44.623\n",
      "Sample:                    01-01-1970   HQIC                           -47.563\n",
      "                         - 01-01-2003                                         \n",
      "Covariance Type:                  opg                                         \n",
      "======================================================================================\n",
      "                         coef    std err          z      P>|z|      [0.025      0.975]\n",
      "--------------------------------------------------------------------------------------\n",
      "sigma2.measurement     0.0010      0.003      0.346      0.730      -0.005       0.007\n",
      "sigma2.level           0.0074      0.005      1.564      0.118      -0.002       0.017\n",
      "sigma2.trend        2.419e-11      0.000   1.61e-07      1.000      -0.000       0.000\n",
      "===================================================================================\n",
      "Ljung-Box (Q):                         nan   Jarque-Bera (JB):                 0.68\n",
      "Prob(Q):                               nan   Prob(JB):                         0.71\n",
      "Heteroskedasticity (H):               0.75   Skew:                            -0.02\n",
      "Prob(H) (two-sided):                  0.64   Kurtosis:                         2.29\n",
      "===================================================================================\n",
      "\n",
      "Warnings:\n",
      "[1] Covariance matrix calculated using the outer product of gradients (complex-step).\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\statsmodels\\base\\model.py:512: ConvergenceWarning: Maximum Likelihood optimization failed to converge. Check mle_retvals\n",
      "  \"Check mle_retvals\", ConvergenceWarning)\n"
     ]
    }
   ],
   "source": [
    "# Load Dataset\n",
    "df.index = pd.date_range(start='%d-01-01' % df.date[0], end='%d-01-01' % df.iloc[-1, 0], freq='AS')\n",
    "\n",
    "# Log transform\n",
    "df['lff'] = np.log(df['ff'])\n",
    "\n",
    "# Setup the model\n",
    "mod = LocalLinearTrend(df['lff'])\n",
    "\n",
    "# Fit it using MLE (recall that we are fitting the three variance parameters)\n",
    "res = mod.fit(disp=False)\n",
    "print(res.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can do post-estimation prediction and forecasting. Notice that the end period can be specified as a date."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Perform prediction and forecasting\n",
    "predict = res.get_prediction()\n",
    "forecast = res.get_forecast('2014')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAD6CAYAAABj5AmjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4XNW19/Hvnq5qFcvdcgf3KlywjQUGAhhCAqaXOJeE3EAIpN1LAn5DCISUG0IJSXAAQwidmG6qQe5N7sa9d1u25aI25Zz9/nFmVEYjaSSNpJG0Ps8zqB2d2bKF9dPea++ltNYIIYQQQojGs7X0AIQQQggh2goJVkIIIYQQMSLBSgghhBAiRiRYCSGEEELEiAQrIYQQQogYkWAlhBBCCBEjEqyEEEIIIWKkzmCllEpXSs1VSuUrpZ6t5brnlVJLlVIPxnaIQgghhBCtQzQzVrcBr2itc4AUpVRO+AVKqWsAu9Z6AtBXKTUgxuMUQgghhIh7jiiuOQEMVUqlAT2B/RGuyQXeDL7+GTAJ2F75AqXUncCdAElJSWMGDhzYwCELIYQQQjSfVatWHddaZ0VzbTTBahEwDfgxsBk4GeGaJOBg8PWTwOjwC7TWs4BZADk5OTo/Pz+a8QkhhBBCtCil1N5or41mKfDXwH9rrR8GtgDfjXBNEZAQfD05yvsKIYQQQrQp0QSgdGCYUsoOjAMidW1ehbX8BzAC2BOT0QkhhBBCtCLRBKvHsJbwTgMZwHKl1CNh17wL3KaUehy4HvgopqMUQgghhGgF6qyx0lqvAIaEvXt92DVnlFK5wCXAH7XWp2M2QiGEEEKIViKa4vWoaK0LqdgZKIQQQgjR7kiRuRBCCCFEjEiwEkIIIYSIEQlWQgghhBAxIsFKCCGEECJGJFgJIYQQQsSIBCshhBBCiBiRYCWEEEIIESMSrIQQQgghYkSClRBCCCFEjEiwEkIIIYSIEQlWQgghhBAxIsFKCCGEECJGJFgJIYQQQsSIBCshhBBCiBiRYCWEEEIIESMSrIQQQgghYkSClRBCCCFEjDjqukAp9UPghuCbacByrfUPwq5xALuCD4B7tNYbYjlQIYQQQoh4V2ew0lr/Hfg7gFLqaeClCJcNB17TWv9vbIcnhBBCCNF6RL0UqJTqDnTWWudH+PB44Eql1Aql1PPBGSwhhBBCiHalPjVWdxOcuYpgJXCx1nos4ASuaOzAhBBCCCFam6iClVLKBlwI5NVwyXqt9eHg6/nAgAj3uFMpla+Uyi8oKGjIWIUQQggh4lq0M1aTsYrWdQ0ff1kpNUIpZQe+BawLv0BrPUtrnaO1zsnKymrgcIUQQggh4le0weobwAIApdRgpdQjYR9/GHgZWAss1Vp/EbshCiGEEEK0DlEVmWutf1Xp9U3Ag2Ef34i1M1AIIYQQot2SA0IbQGuNYda0KiqEEEKI9kqORWgAU8PxIi9Ouw23w3o47JJRhRBCiPZOglUj+A0Tv2FS5AW7TQVDlh2XQ0KWEEII0R5JsIoRw9SU+AxKfAZKgdthL5/NUkq19PCEEEII0QwkWDUBraHMb1DmN1CA22knyWWX5UIhhBCijZNg1cQ0FSHL47ST7HZgt8kMlhBCCNEWSbBqRmV+A6/fwOOyk+xyYJOAJYQQQrQpEqyamQZKfQZlPoMEl50kCVhCCCFEmyHBqoVooMRnUOo3SHQ5SHLZpchdCCGEaOUkWLUwraHYG6DEFyDJ5SBRApYQQgjRask2tTihNRR5Axwv8hEwzJYejhBCCCEaQIJVnDG15nSpH62lZY4QQgjR2kiwikMBU3PWG2jpYQghhBCiniRYxalSn3X2lRBCCCFaDwlWcexMqV/qrYQQQohWRHYF1pPfMJm1YBe7CorITHaTkeQiM8lFZrKLzCTr7Vg1YdbA6VI/GUku2SkohBBCtAISrOrh4KlSfvTqatbsO0Wqx8GZssh1UKkeBxlJLit0Jbvp2sHDlcO70iszqd7PGTA1Z8oCdEhwNnb4QgghhGhiEqyi9NWWY/zkzbUEDM3TN41ibJ8M/IbJyWIfJ4t9nCjycaLYG3zp40SRlxPFPjYfPsNXW47xr6V7mdg/k5vHZjOmV3q9ZqDK/AZuhw2P096EX6EQQgghGkuCFbB06VLy8vLIzc1lwoQJVT4WMEz+/Pk2/p63k0FdU/nbLaM5uHUdT/35ec6fPJmcsePpnOqp9f4nirzMWX2Q/6w+wN2vrmFAp2RuGpvNJYM7R71seKbUj8OmcNilLE4IIYSIV6olzkvKycnR+fn5zf68kSxdupSpU6fi8/lwuVzMmzevPFwdOV3Gj19bw4o9J7lpbDa/vmowa/JXlF/vdLl4+/255IwdH9VzeQMGn248ymsr9rHreDGZSS6mj+nBNaO7k5boqvPzHTYl9VZCCCFEM1NKrdJa50RzbZ3TH0qpHyql8oKPtUqpZ2u47nml1FKl1IP1HXBLysvLw+fzYRgGPp+PvLw8ABZuL2DaUwvZeOg0T9wwkseuGYbHaa9yvd/nY8nChVE/l9th55sju/Hq98fx5I0jOadzCs8u2MU3/7qYx+ZuZvfx4lo/P1RvJYQQQoj4VOdSoNb678DfAZRSTwMvhV+jlLoGsGutJyilXlBKDdBab4/5aJtAbm4uLperfMbqgilTePzzbTz95XYGdErmb7eMoX+n5IjXO10uzp88udb7569YxpKFC8uXDQGUUozvm8n4vpnsKiji9ZX7mbvhCO+uPcSEfpn8cEo/zu2SEvF+Um8lhBBCxK+olwKVUt2Bv2itr4/wsaeAT7TWc5VSNwIJWuvZNd0rnpYCoaLGauT4C3h5h50lO08wfUwPHr56CImu6tlz0eIlfPzZvCphKZL8FcuY/s0r8EexbHiy2Mec1Qd4e9UBtIZ/3TG2xtotBWQkuaTeSgghhGgG9VkKrE/x+t0EZ64iSAIOBl8/CYyOMKg7gTsBsrOz6/G0DVNbQXq4CRMmoDsN4N7X13K2zM8fpw/n+pyetV4/YFi1L7GaJQsX4g8uGxJcNqwpWGUkufje5L5cMrgzM2av5JdzNvCPW8dELG7XwKlSP5n1qLfSWhMwNXalsNmkRksIIYRoClFNeSilbMCFQF4NlxQBCcHXkyPdV2s9S2udo7XOycrKasBQoxcqSJ85cyZTp05l6dKltV6/YvdJbn1uOSkeB+/ePbHWUFUf50+ejNPlwm63R7VsCNArM4mZVw7m60NneHJezaupRg31VlprfAGTUp/B2TI/p0p8FJz1cuysl5PFPo4Xeyn1SascIYQQoilEO2M1GViua143XAVMApYBI4CtMRhbg0UqSK9t1urrQ6cxNbz2/bqPTqiPnLHjefv9udVqrOpy0cBO3Dwum1eX72NY9w5cNrRLxOvK/AY2Zc1gGYY1I2XWsbSrNZwp81PqN0jxOHDKcqIQQggRM9EGq28ACwCUUoOBm7XWlXf/vQssVEp1Ay4HoksQTSS8ID03N7fW670Bqx9fiif2x3rljB0fdaCq7O7cfmw6dIbHPt7MgE7J9KtUQF9ZSQNnn0KHm3qcdlLcDlkeFEIIIWIgqukKrfWvtNZzgq9vCgtVaK3PALlYM1YXaq1Px3qg9TFhwgTmzZvHb3/72yrnUtXE67eClSuOZm8cdhuPfnsoSS4H98/ZQJG3aY5ZKPMbHC/2UuKTYxyEEEKIxopZktBaF2qt39RaH4nVPRtjwoQJ/PKXv6wzVAH4DCMuTzXvmOzm0W8P5WBhKb/9cBNNdZir1nC2LMCJIi++4OydEEIIIeovvpJEC/H6TdxRtpZpbqOy07n7on7kbS3g1RX7mvS5AqamsMTH6RI/ptn8J/ILIYQQrZ30CsSqsXI38sBNpcCuFHabdZyBw6awBd9Wweco9RsYDQgsN4/NZsOB0zzz5U4Gd01lVHZ6o8Zal7KAgbfIIMntIMkt3yJCCCFEtOJzmqaZeQNGveqrFFahe4cEJxlJLrKS3XRK8ZCZ7CYt0UWqx0miy4HHacdpt+Gw20hyO+iY7CY90UWCy0592v0ppXjwysF0T0/ggXc2crzIC1gHkD715z+Rv2JZPb/iummgyBvgZLGPgCHLg0IIIUQ0ZDqC0IxV9MHKZlMRT2SPhsthw+WwkeJ24A2YlPkNfAGTuuaxkt0Ofn/NMP7rpZU88M5G/muAjxu+NS2qU90bI7R7MNHtIFlmr4QQQohayYwV4As0f42VUgqP005aoouOyW5SPA4cdRx50K9TMr+8fBBr95/i7/N3l5/qXt9m0PWlgWKvVdzul9krIYQQokYyBUFwxsrRck2NQzNgiS4HCxct5tMvvmTCpMgHil42tAsbDp7m7VWQPPgCijYtqPVU99OlfrYfPcv2Y0VsP1ZEtw4ebh6X3aAZt4Cprdkrl51ktyPqdjpCCCFEeyHBCqvGKh52BS5dupRvXHoJvjqW9+67eACbD59h51U/45KrL+OyKWMZfd449p8sYfuxIrYdPcuO4MujZ7zln5ee6KSwxM87aw5y94X9uWxoF2wNCEclPgNvwCTV44zYy1AIIYRoryRYYR23EA8BoXIrntqaNjvtNh67Zhi3P7+CVe4x7N3iYMeC+eWnsNuU1XNwRI80zumcwoDOyQzolExmspsNB07z+Ofb+M0Hm3h71QF+esk5DO3eod5jNYJHMyS4rJPbZfZK1EVrTbHPoNRnkOiyk+iyy/eNEKLNkWCFtRTYFO1s6iu8Fc+FF+bWeG3nVA+PfGsoD3+4CQVMG9aVAZ2TOadzCn06JuGp4fiIYT068PyMHD7ZeIS/frmDO17K5/KhXbj7wv5kpbjrPeZSn1V8n+JxtOhyqohvZX6Ds2WB8l6WRd4Axb6AtQTutEtLJSFEm9HyaSIO+Fq4xiok1IonLy+P3NxcJkyYQGGxD18NBePn9cngg3smRX3//BXLyhtCXzF2PFPOyeKlpXt4dfk+8rYWMOP83tw0rme9/ywMU3OqxE+iyyTF46zX54q2zW+YFJUFIn4Pa21tiijxBkhw2Ul0ObBLwBJCtHISrAjWWNXjuIWmNGHChCpteDokODlR7Cv/Tb+h8lcsY/o3r6h2PMNduf25ekR3npq3nb/P38m7aw9y79QB5J6bVe9lmhKfgd/QdEhwyg/Ids40NUW+AKVRNAnXWN87pT4Dj8tOkgQsIUQrFh9pooV5W+C4hWjZbIq0RCeN/TGzZOHCGo9n6J6ewB+mD+evN40i0WXn/jkbuPvVNWw/drbez+M3TE4Ue/EG6v6BKtqmEl+A48XeqEJVZRprafl4kZfTpX45mFYI0SrFZ5poZt5AfBSv18Rpt5Ga0LgltvMnT8bpcmG322s8nuG8Phn8646x/OIb57L92Fluf34FT36xvf4/IDWcKvFT5A00asyidfEGDE4UeTlbFqCx/cLL/AYnin2cKvGV1/E1VRNyIYSIJVkKBLx+Iy5qrGrjcdrxG2b5zr/6yhk7nrffn1teY1XTKe0Om43pY3pwyeDO/D1vJ6+u2EfetmP86vJBnNcno17PWewN4A+YdEhwSnFyG2aYmqKyAGVNMEvpDZh4AxUzVzZl9eF02BUOmw17qC+nfH8JIeKEBCvAZ8TvUmBlKR4nAUPXWMxel5yx46Nue9Mhwcn9lw/k0sGd+d3Hm/nRa2u4cnhX7p06oF6zZz7D5ESxjw4JcuZVW2OYmiJvgEWLF9cZ2GPF1BqfobF+v6gIcpUDl9Nuw2W3SdgSQrSIdh+sDFPjN3Tcz1iFxKqYPRqhXYQ/PX8Sa8p68cqyfSzZeYKfX3oOFw3sFHVxu6k1p0p8JHscUZ/4bppWgPQZJr6AicOmSPXIzFc8MExNsS9Amc9gZQ2bIppbpMBlt1WELKdd4ahHo3UhhGiodh+sfMFlhnjZFViXUDF7YbGvzsbNjRFpF+HF3z2PR+du5lfvbGTKOVn84hvnRn32lQbOlgXwBzSpCdUPFK0cpPwBk4BZ9aszTC0zXy3MDAaqUp9R/r1XeVNEbYfalt9DawqLfWQm1//MtPoyTI1hGpT5raClFLjsVhN0p916CCFErLX7f1lCu9dcregf2foWs9uUwu2wkex24IlyZi7SLsJzu6TwwowcfnRhf5btOsENs5byzpqD9Zo9KwtYRcm+gEmZ3+BMmZ8TRV4KgjvBSn1GtVAVYmrrtHcpiq8/w9QEjIYVgJum5myZn+NFXkoqhSqIblNEyMaDp/neS/lc8dQifvjvVazYfbJZC9K1tmq2zpYFOFns49iZsvLQJYQQsRL1jJVS6m/Ax1rrDyJ8zAHsCj4A7tFab4jNEJuWt5XNWIV4nHZ8hlltx55S4LTZcDpsOIJLIeFnAiUE7BR5A/hrqdUK/cAkOGMV+oHpsNm4bUIvppybxWNzN/P7j7fw2ddH+NUVg+iZkRjV2EPtcBqq2BvAFyyKj+fzjnwBE1PrGk/Bbw6RZplCBd8Ou/U9Eno7fBYx1IKmxFfzLr9oNkUcL/LyzFc7mLvhCB2TXdw6PptPNx7lntfWMKRbKt85vzeTB3RsUN/KxtBYTcq1hgRX6ygFEELEPxXNb4xKqcnAT7TW19Tw8dHADVrr/43mSXNycnR+fn69BtpU9p8sYfIfv+L/rhvB9DE9Wno49Xa61F++xBH6YRmt8DYj4Sqf1B7pB6bWmvfWHuLpL3fgN0x+MKUvN43NbrYfkEpBqsfZosGlMq11+S42b8AoDyN2myLRZSfB2by98Up8AYq80R99EApY9mDIqi1QRcMXMHlj5X5eWLwbv2Fy09hsZpzfmyS3A1/A5KMNh/nX0j0cOlVGv6wkvnN+b6YO6oTD1vy/5KR6nBKuhBA1Ukqt0lrnRHVtXcFKKeUENgBzgfla6/ciXHMXcDdQHLz2B1rrGtdr4ilY7Th2losfX8DTN43iqhHdWno4zU5rTYnPoNgbaHDNVsFZL3/8dAsLth1nUv+O/L+rBtOhjqXKukJbfXicdlI9LdMI2jA13oB1zpIvYNb6Z6gUzdIbr8xvUOQNYNSwpFqXWPzdLNpxnL98vo0DhaVMHtCRe6cOoGdGYrV7B0yTzzcd5aUle9l9vJge6QncPqEXVwzr2uw1UCn12FwhhGhfYh2s7gCmAXcB9wBHtNZPh11zHnBAa31YKfUv4G2t9fth19wJ3AmQnZ09Zu/evdF+PU1q48HTXPn0ImbdNoZLh3Rp6eG0GNPUnPUGGlxzorXmrfwDPDlvOx2T3Tz67aEM7d4h4rU1tddpDLtN0SHB2Sw/jP1GcFbKX3M9WG0U4HHZSXTaY7pTrba+fNFq7N/N3hPF/OWL7SzdeYJeGYn89NJzGN83s857m1qzYFsBsxfvYcuRs3RKcXPLuGy+Nap7s85ISrgSQkRSn2AVzb8go4BZWusjSql/A48CT4dds15r7Q2+ng8MCL+J1noWMAusGatoBtccQjVW7X2nmS0YTBJd9gb9cFZKcf15PRnSPZUH3tnID15exY+nDuD6nB7VZpLqu5MsGoZp7Tarz5EOdd2v/KE1hmG9DJhmo08VD7VuKfUZuB02El2ORn3/NTYUV9bQv5sib4AXFu3mjZX7cTtt3Dt1ANfl9KgSdGu7t00pcs/txJRzsli++yQvLt7DX77Yzj/nb6erUcDIAd3JHTOEPh2TyEhyNdnsZOjU+CS3hCshRMNE86/HDqBv8PUcINJU08tKqUeBjcC3gN/FZnhNL7QrsLWcY9XUnHYb6UmuBi8nDenWgZf+ayy//XATj3++jTX7Cnlw2mCSPRXfajUVxjdW6EiHMr8ZrBOy3h/6ERz6YVzxdvDzNFXCk2nqJj3KojKrHsuH024j0WUvLyKvPMZIIWLp0qV89dVXjJ84maGjzovZeOv7dxMwTeZuOMLf83ZSWOzjyhFd+eGUfhGPU4jm3kopxvfNZHzfTF7/bDG/f3MRZ7oNZPtWH29tXQNAqsdB745J9Al7dEpxxyRwFQWXxZMlXAkhGiCapcAU4AWgM+AEfg5crrV+sNI1Q4FXsX5mva+1fqC2e8ZTjVXe1mPMmL2SOXedz+js9JYeTlwp3xUWRf1VeO2M1ppXV+zjmS930qWDh99dM5SBXVJrvL61aKlxq+B/FIr8Fcu45qrLm+xQzmi+RlNrPt90lOcW7mbfyRKGde/Azy49h0FdUyNeX597hzz15z/xh0d/g2EYOFOzuPXHv2Jk7jR2Hy9mz/Fidh0v5nSpv/x6jx2SzSKG90zj0pxBDO/RoVHnZSW5HRKuhBBAjJcCtdZngevC3r047JqNwPCoRxhHyo9baOdLgZEopYJnX9k4U1bz8Qw11c7cMq4XQ7t14IF3N/L9l1bxk0sG8O1R3VFK1dleJ2CabD58llV7Csnfe5Ikt4OfXnIOnVM9TfXl1qkpasOipYP/0WgWLVgQ86XUymr7u9Fak7e1gFkLdrHreDH9s5L54/ThXDCgY1SzRfVpq1R5hsvuL+LaKSPJOa9nlWsKi33sPl5M3qqNzH7jPU5nZFPgd/HlPuu0l54ZCYzokcaInmmM7JFGz4yEqGe1ir0BtNakeBrXAF0I0b60+1/HKoKVLAXWxGG3kZHkotRncNbrr1ZjVFvtzIieafz7jrE89P4m/vDJVtbsO8X9lw+sVsNias2OY0XkB4PUmn2nyhtO9++UzMaDZ7j1+eU8eMVgppyb1Sxfd7imqA1riIYspTZ2pk1rzeKdJ5i1YBdbj5ylV0Yij3xrKFMHdWqy4zWiOSMrPclFepKLRW8v4+Tn/8AwDOxON//1v7/jnMnTWLf/FAu2F/Dh+sPW9YlORvRIY2R2GiN6pHFOl+Raj3cIfQ9KuBJCREuClT9UYyUzVnVJcNlxO2xWHVOgolC6rh/0aYkuHr9hBP9aspdnF+xk65Gz/O6aYTjtKhikClm1t7B8WSc7I5HLhnQhp3c6Y3qlk5boYt+JEma+t5H/+c96vj2qO/ddPCAmu8XqEziaqjasvqIJHJWtWL6UW++diWHCk/98iVdf/Cfjxk+I6rm01qzcU8izC3ay8eAZuqV5+PVVg7l0SOdmOW8q2hmuKn83DhtXTxlNzthe3Da+F6bW7D1Rwrr9p1h34BTr9p8mb1sBAIO6pvDkDaPokFhzcAqdNp8q4UoIEYWoDgiNtXiqsfr3sr08+O5GVvxqKp1acJmptQk/XDTagLJqbyEz393IieKKk9c7p7rJ6Z1BTi8rSNW03Oc3TP4xfyf/XraPPh2T+O23hjCgU0qDv4aGLO1F+3WW+AIs33WSU6V+LhnUuUrxfnNat/8Uv/z3Ak7opPL32TDp3TGlvOg7VAienZFYZXfimn2FzFqwi9X7TtEpxc0dk/pw5fCujToiwuOw43HZKPOZVcJ5LNQnJBec9bJ4x3H+/Nk2sjMTefqmUWQkuWr9HHewx6C9ltPqhRBtU0zPsWoK8RSsXli0m4c/3MS6/3dprb+1iuq0trb5h7fVqcuJIi9v5O+na4cEcnql0yO99rqX8B+Yy3ef4Dfvb+JsWYAfXdS/2pEO0f6ArVwcbbfb+d8Hfs2Pf/aLen0tlRWc9bJwewELtx8nf09h+ZEViS473x7VnRvO69lsNWKHTpXyzFc7+GLzMTq4Ffs+eJqyY7vxdOrN9O//lDJXB3YfL+bQqVJCGz9tCrqnJ9CnYxKlPoOVewrJTHIx4/zeXD2qW6OWyx02RYqnagPtgGFS7DPw+o1m24UZbsXuk/z8rXV0S0vgmZtH1bvYXSmrzZNdKex2Zb20VTyEEG1DrM+xatNaa6/AeKCUslqBOO2cKfVHfVhmZrKbu3L7R3VtpFmlcWPH88r3xvHIR5t5/PNtLNt1gplXDiYjyVWvWajGLu1prdlVUMyC7QUs2HacTYfPANAtzcM1o7vT2Shgy4a1nEgfwusr9vP6yv1cMrgzt4zL5pzO1WfaYrHjsMgb4KUle3h9xX6Ugu9N6sOt43vx9eTEiPcu8xvsO1nCnuPF1m67EyXsPl5MiS/APRf1Z/qYHo1aclUKUtyR28U47DY6JNgw3A6KfQHKfM0fsMb2yeAvN4zkZ2+u47//vZpnbhlFp5Tow6/W1kyqHyCs14Tdpkhw2vE47RKyhGhHJFgFlyNczdw+oy1x2m1kJrsp8gYo9tbYyahBaioYT09y8X/XDeftVQd4at4ObnluOb++ajAr6lFgXt9aJbB2K67bf5oF2wpYsL2AQ6fKABjSLZX/ntKXCwZk0TcriVUrlzP9m1eVB7xn3/iQzf6OvLf2EJ9sPMLY3hncMj6bcX0yUEo1esehYWreX3eIZ+fvpLDEzxXDuvDD3H7lIaGmWiWP0845nVMiBr3GSnDZSXY56mzfY7dZAT3ZZQWsUr/R6ENY62NMr3SevHEk972xlh/+ezXP3DyaLh0aP7NomJoir9Wv0WW3ldcoyvKhEG2bBKuAidOumrR3W3uR7HbgtCtOl1bfOdhQtc0qKaW4LqcnI7PTmPnu19z7+lqmZo/F6UmEspKoZqEiBY6AaXLsjJcDhaUcKCzh4KlSDhSWcrDQelnqN3DZbeT0Tuf2Cb2Z1L8jWSlVl5DCA+HW/MXc97NfcMekPry75hBvrNzPva+vpV9WEreM68W2BYsavONw+e4TPPXFDnYUFDGiRwf+ckPd50k1JZfdRorHUe9aLFtwuTDJ5aDUb1DcyCbQ9TGiZxpP3zSKe19fy3//exV/u2U03dISYnZ/n2HiKzVRygqzCU57s/dCFEI0j3ZfY/XwB5t4M38/G3/zjZYeSpthmJozpf5G9ayrLJolsjK/wVPztvOf1QfpkWwjpfQwPXr0oEf37uX1LrZKRcd2VfF2qc+oEqAOny6rcuK8067onpZA9/QEuqclMDo7nXF9M2ptnVPXDJTfsJoPv7JsHzsKiujgVhz48t+UHdiMXWn++H9/ZtTIkbgcNtwOG26HHbfTVqVges/xYp76cjuLd5ygW5qHey4awIXnZrXYjIhNKVI8jpj19tNaU+o38PpNAqYu3yjRlDYfPsOPX1uDx2nnb7eMpmdGYpM9l8OmSHDZ8Tiatim3EKLxpHi9Hh58dwMfbzjCqpmXtPRQ2pymWBqsy/xtBTz+2TZOlfowTI1pWu1q6pLsdtD3G9EWAAAgAElEQVQ9PYGe6VaA6pGWaL1MTyArxd2gs5qiCYRaa5bvPskry/exYvfJOu+psOoBXQ4bxWUGHpeN707sww05PVus36UCEt0Oklz2Jg11Wmv8htW/MWCaGKb1dqwD17ajZ7nn1TU47Ipnbh5N745JdX9SIyisc/QSXPZ237NUtDFaWw8ivAx9vNb3UcPHqXptQgY08fEvEqzq4RdvrWPxjuMs+eXUlh5Km+QNGDFdGmwIrYONlENBq3JzZVPjcthI9ThavPZl/8kSjhd5g/0DTbx+A59h4vWb5S+9AQNvwMQXMEl2O7j+vJ51HhPQVBTgibKOqilprQkE/079honf0AQMs1GF8LsKirj71TVorXnm5tH065Qcs/HWJjSLleCMTUg1TY3PMMuPiRCiXHjo0SYRA1D5S7N66KkzNDWT5M5xFaykxipg4o7R0oWozu2w0zHJxulGLA1WbpbcsM9XOJQi3icDemYkNunSU6yEAlWSyxEXP6yVUjjtCqed8mXI0OyWFbSsYBr+/VPbjGLfrGT+fsto7n51NXe9spqnbx7VJAX+4QKm5myZVfDucdpJdNrrXatmmJoyvxXAK7ehctptJDitAvrGBmGttRXwDROnzRaTe4YEDJOygElZ8PBmR3AZ31G+lN/OQ6IOD0JmhLcrvw9qDk6iKUiwChiyI7CJ2WyK9CRX1EuDCnA5rOUul92Gw26Led2WaBiP006yOz4CVW2UUrgcqtq5WX7DmsFZsmRJnbswe3dM4h+3juHuV1dz9yur+eFwBwfWLWmWBtxaQ6nPoNRnlO8orK12zW9UzHLWdOxJKGSG/v/yOOu3SzEQfA5fMLCFnqWUip3VngYGN9PUlAWsrzd8/EaEr0dBedgKvWyVtWqmWSkIRXhUnk2qHJhEXJNgFTDlDKtmUtOuQYV1plEoSEWqM7EHw1lxMJzJPy3Ny+Owk+Su/+xJPHHYbTjskICdtSuWRLULs2dGIv+4dQzfe2Epjy06g/eAyT+Xv8Xt1/sZP2Ig3dMS6JaWELOC/UhCOwrPlgVIDC4T2mwKX8BaGi7zm/WqM9NQvtysVLC+y1m9vktrK4SGwlSkgFNtnMFffKIJWaFZr1KfUe9fmDQVS/oh9kSF2xYHqw9ag2kEQ1DwZfnb4Q/5l6wtavfByhcwpU9gMwotDZ4tC2CzWcsT9fmtOcntwOWwlhbr+odeWNwOG8luBzalymuRoi3+djtsJLkdbe5ogNzcXFwuFz6fD1cdx3J0S0tgsrmeV3eexJHRE0ePIby9zcfb29aXX5OV7C7fNRp6ObhbKtkxXNo1ta6Y9VWx+ZmstbWjtsxvYFMKj9NaZvMFw1RDn6K2kFUW3OnpDbTcifuNYgYDk2mAGah4vXKAEu1auw9W3oBJgtRYNSubTTWqfZDTbiMzycWZskB5HUY8CYVFt8OGTSlK/AYlzXgmU4jDpkj2OKq0onGVzx5UvK9y8XfAtAq/gTYZqEImTJjAvHnzyMvLIzc3l8EjcmrtXXjJBeN59v8qlg5nv/0RWX0Gc/BU8HyzU6UcKixlxZ6TFGzwln/e8B4duHJ4V6YO6kyyOzb/3Ory/8SWqTUl9WxPFY3KIUvFKBA2Oa2t0GQGwPAHXzeCM1Ct4QsQLand7wq88umFdE7x8PyM81p6KKIByvwGZ8padtdhqGbF7ah56SN0JlOJz2jymbba2siImp0u9dca1KNtOeQNGHy2YDkfrdzOYU9PjhRrPE4bFw3sxJXDuzEqO61Bx3eI2qUlOhvWz9I0KoUnPxjBQCVaD9kVGF+8flPOjmnFPMETrE+X+qvsgGpqSoHbbh3aGc1SplKKRJeDRJeDsmDAivV4FcE2Mu6WPzqiNeqQYM2i1hSuamoLFG7D6pXcc/OV5bNbf3r5A3b4M/h881HmbjhC97QEpg3vyrRhXRvUOicWPSXb2r3rzTSg9FRwKU+W7kRsSbCSGqtWz25TZAQL24ua4EBSFXyO0CM0O9VQnmBjXl+wcLe2Jaio7+mwk+yJ/9168a5DghOlrB15DRXezujw+iX88me/4CeXnMNXW4/x4brDzFqwi38u2MV5fTK4anhXppybFdX3VGN7SrbFezeINsHwtdzzizat3Qcrq3hdlkzaglBN0Jmy+he2KwV2FTwjx261vLFX2s7dFEJHSiSbVvPhMl/9i3kdwf56MusaO6keJwoaXG9UU39Lj9PO5UO7cvnQrhw6VcpH6w/z4frDzHzva1I8DqYN68q1Y3rUWvBeU1PyWGit9xYi3kQdrJRSfwM+1lp/UMPHnwcGAx9prR+J0fianDdgyHELbYjLYaNjstUQObx+MLwOK/Rm6EyclmK3KVI9TlLcDgxTo6l0sDG6yiHHOjhqra3Pa8pt/u1ZisdaFmxIuMoZO563359b67JXt7QEvn9BX+6Y3IdVewp5b90h3lp1gNdX7md83wymj+nB+f06VpuBrK0pebjTJX4W7TjOqr2FDOqawsWDOpNeyyn99bl3fdX33lprPvhqKetX5XPJ5PMYN35CzMYiRFOLqnhdKTUZ+InW+poaPn4N8E2t9Qyl1AvAY1rr7TXdL56K1wf/v0+4ZVw2D0wb3NJDEULEmebqd5m/YhnzFizF12MMy4/ZKCjy0rWDh2tH9+CbI7pV2UVbW63S0TNlzN9aQN62AtbuO4WhNcluB0XeAHalGNc3g8uGduGCAVkRNzdEUwd1ptTPsl0nWLzzBGv2FdKnYxJTB3ZmyjlZte72jebeBwtL+eTrI7yXv5ujJcFfIkyDrCQH/bum0zM9kR7pCfTISCQ7PZGuaZ7ynav1Kl43/FB8PLprRfyLs+L1OoOVUsoJbADmAvO11u9FuOYp4BOt9Vyl1I1AgtZ6dk33jKdg1e9Xc/nvKX35xTcGtvRQhBBxqKnDVXj90RvvfkRRh768veoAq/edwmW3ccngzlyX04NBXVOrfK7Wmt3Hi5m/rYC8rQVsOXIWgN6ZieSe24ncc7Mo2vc1cxfk4+82gvWFNo6e8ZLgtDPl3CwuH9qFnN7pOGr5oaS1ZldBMYt3HmfxjhNsOHAaQ2s6JDjpk2yyo6CEItzYbYrzeqdHFbIqO1Xi44vNx/j06yOsP3AagM7qLFs+/ReG34srowfDJl5MQqde7C8sodhbMYtoU9Clg4ee6Yn065TMg9MGkeiKYiFGglXbEmfBKpqlwNuBTcAfgXuUUtla66fDrkkCDgZfPwmMjjCoO4E7AbKzs6MZW5MLGNYhiVJjJYSoSbLbgYIm2RgB1euPli9exI9/NoGpgzqz81gR/1l9gLkbjvDRhsMM6ZbK9DE96JmeyILtVpjad7IEgCHdUrn7wn5MOSeLXplJgBXarrt6Wnloe/O9j3B0Gc0nXx9h3uZjfLLxCOmJTi4d0oXLhnRhUNcUlLIO8czfW8iSHcdZsvMEh0+XATCgUzK3TejFpP4dKT24ieuD907oPpAb73+cr0+W8Ojczfz+ky3k9Ern4kGRQ1aZ32DR9uN88vURluw8gWFq+nZM4u4L+3Hp4C4c2LqW6U/mWX8uLhf/M/O75Iwdi9aaUyV+DhSWsr+whP0nS8pf/+zrIzxy9dAm+TsSoj6iCVajgFla6yNKqX8DjwLhwaoISAi+ngxUi45a61nALLBmrBo84hgKHVonuwKFELVJCvZH9BkmZvAgVdPUMTmns7b6o36dkvmfywZyV25/5m44zNurDvCbDzYBVo3dmF7p3HheTy44J4usFHe1e4eHtmWLrNA2ulc6P7v0HJbsOMEnXx9hzuoDvLFyP9kZiXRPS2D1vkK8AROP08Z5vTP4zvm9Ob9fJp1TK46HeOqNReX3Lj24hU5HlvPbn/6cLUfO8uWWY8zbfMwKWR9vIad3OlMHdaJzqofPNh3lqy3HKPEZZCW7uWlsT74xpAsDOiWXHxPSpYY6NaWs1lbpSS6G9ehQPpb8FctYvXwFy5cnM2GC1GOJlhVNsNoB9A2+ngPsjXDNKmASsAwYAWyNyeiamNcvwUoIEZ3QMRmVhfrVGabG0BrDsF4GTDPqQ2ujKXZP9ji4/ryeXJfTg1V7Cyks8TOuTwapCbUvt9UW2twOOxcO7MSFAztxptTPV1utGaxDp0q5emQ3JvbvyKjstBpn9CPdWynFoK6pDOqayl25/dh69CzzNlsh63dztwCQ6LJz0cBOXD60C6Oy02s8IiTac8MqL6U+9uijzJs3T8KVaFHRBKvngReCtVNO4OdKqUe01g9WuuZdYKFSqhtwOdAq9tF6A8FgJTurhBANEDrbLJKAYVLqNyj1G3WGrGhDhFKKnN4ZUY8vmtAGkJrg5OqR3bl6ZPeY3VspxcAuqQzsUhGyjp/1kdM7Paa7WSvPyvl8PvLy8uoOVtKWRjShOoOV1voscF3YuxeHXXNGKZULXAL8UWt9OmYjbELe4MGMrjbaD00I0XIcdhspdhspHmd5o+PGNDZuqGhDW1PeOxSy6BL7MVSeOXO5XOTm5tb9SX36QXIi9OsD/fpC/74VLzt3sg62E6KBYnZAqNa6EHgzVvdrDhUzVhKshBBNJ7SMaJqasoBBqc8g0MQ9I9uL0MzZ6uWLufiii+qerTIMuOF62LoJduyCz78CX6VT2JOToG8fK2RVDlx9+0BGetN+MaJNaNcnr/tCwUp2BQohmoHNVtEz0h9cKiyLYqlQ1C5n7Hguzp0c3b/ldjv84fcVxy0YBhw4BDt2wq7dVtjasQvWrId3PwSzUi/BjPSKsFU+0xWc9Uqq+cR80b6062AVWgqU4nUhRHNz2m047TZS3A68AbNJGnOLKNjt0Kun9ZiaW/VjPh/s2WsFrZ27K17OXwSvvVX12m5dw5YWg7NevXuBq+YT70Xb076DlewKFEK0MKVU+VJh+SxWA/pGRmJTCpfdht80690/U2AFonMGWI9wxSUVM1w7d8POXdas1/sfwcnCiutsNsjuWTVshWa8ena3gp1oU9p3sAouBUoDWyFEPAjNYiW7HJQFDEp8Rr0Dkd2mcDtsuB32Kv+2+QImZYHYLD0qZZVQuB02/IZJmd/EjNF6ptNuw+O04Qku6wVMjamts8MacpxFk0lKhGFDrEe4k4VW6Nq52wpbofC1bAUUFVdc53JB395VZ7r697PeliL6VqudB6vQUqD8xiCEiB+Va7G8wWL30C+CkTiCDbndDhuOGnY5uxw2XA4bqR4n3oBBmc/EG4h+ZsxhU7iDz+Gs9Bwep50UD8GAZTQoZIUaiic47dWOr3DVcJyFGTo7LHhgqz2eQkhGuvXICWtCojUcPRYWuIKhK1IRfb9Ks1z9+1UU0aenNe/XI+qlnQcr2RUohIhv1syQHcPUlPgC5ediuew23E5rZqqms7TquqfWVn1Xqc8o70QRorDCWGhmylbHc5TXjEUZspSiPEw5G3Dkjc2msKFoVccQKgVdOluPiWHHVFQuot+5q2KWa/W66kX0mRmRi+j79pEi+jggwQqpsRJCxD+7TZHicZLsdqA1dQadaFSu7wodBWGY2prdstvKW8zUV00hS2urN6vHZZOVgnC1FdF7vbB3H2zfWVFEv2s3fLUAXg075ah7V+gbVkDfvy/0ypYi+mYiwQpZChRCtB5KqSYpvQktP8Za5ZCltW5wWGvX3O6ai+iLiiMU0e+C9z6sWkRvt9dcRN+jmxTRx1D7Dlb+4MnrMmMlhBBNTkJVE0hOguFDrUe4kyeDYSt0PldwxmvJMmtXY4jbXb2IPvSyU5YU0ddT+w5WshQohBDtj7KDMwECZW27b2BGhvU4b0zV92sNR45WOpsrGLq274TPvqxaRJ+SXL2IPvR2Wofm/XpaiXYdrHwSrIQQov2x2SAhDUwDfMXgL2nbASucUtC1i/WYFNYCyDBg/8FgHdeuipPoV62Fdz6oWkTfMbPS6fP9KhXR94bE9ltE366DlTdg4nI0vEBTCCFEK2azgycVXMngLwZfCeh2fvq93Q69s60HuVU/5vXCnn3WDFeokH7XbvhyAbwSoYi+ctiqXETvdDbXV9Mi2nmwMmS2Sggh2jubDdwpVsAKzWCZRkuPKv643XDuAOsRrqi44kyuysdFzHkfTp2quM5ut8JV5T6L/SsV0dta/8/kdh6sTAlWQgghLEqBOxlcSeAvtUKWGWjpUbUOyUkwYpj1CHfiRFgRfbCua9ESKCmtuM7jgT69wvotBme9sjq2miL6dh2sfAFTjloQQghRlVLgSrQeoYBl+Ft6VK1XZqb1GJtT9f2hIvrKYWvnLti2HT79AvyV/sxDRfSVA1eoiD65c/N+PXVo18FKZqyEEELUypkQ3EHos5YI2/pOwuZUuYh+8vlVPxYIWEX0obAVCl8rV8N/3qv6dzByJKxZ07xjr0X7DlZ+Q86wEkIIUTeHy3qYphWwpA6raTkc1rJgn15w8YVVP1ZWVlFEv2MXJMTXsQ/tO1gFTNytqtGUEEKIFmWzWXVY7mTwlwVnsbwtPar2xeOBgedYD2hdS4FKKQewK/gAuEdrvSHCdWuBUNn/o1rrz2M6yibiDRi4G9D8UwghhMDpsR5GIDiLVSrHNYg6Z6yGA69prf+3pguUUpnAFq31jTEdWTPwBUyS3O160k4IIURj2R1gT7WObPCXWg/DV/fniTapruma8cCVSqkVSqnngzNY4cYBY5VSS5RS7yqlUmI/zKYhxetCCCFiJrSbMCkTkrKsYxuU/Ixpb+r6G18JXKy1Hgs4gSsiXLML+IbW+nxgPfDdSDdSSt2plMpXSuUXFBQ0Zswx45XjFoQQQjQFu8M61T25EySkg8Pd0iMSzaSuYLVea304+Ho+EOG4VXYBO+q4Bq31LK11jtY6Jysrq0GDjTU5eV0IIUSTUsqqw0rMsEKWOwVsUoLSltWVKl5WSo1QStmBbwHrIlzzKHBV8PXpNVwTl7x+U45bEEII0Txsdms3YXIWJGZagauVnCYuoldXqngYeBlYCywFViulngu75nHgAaXURsALvBTzUTYRnyE1VkIIIVqAw2UtESZ3Bk8HsLtaekQiRmqdj9Rab8TaGVjZ98KuOYxVwN7qeP1yjpUQQogWVLl9jmkEj20okx6FrVi7XejVWkuNlRBCiPhhs1s1WO4Uq4VOoNQKWXI2VqvSboNVwNSYGglWQggh4k+ohY471TrZ3V9inY0lfQrjXrsNVt6A9RuAFK8LIYSIW6FdhU6P1acwNIslB5DGrXYbrHzBYCXnWAkhhGgVbDbr0FFXUrAeK3jKu9RjxZV2G6y8AasruSwFCiGEaHVCRze4k8HwWwErUGYFLtGi2m+w8gdnrJwSrIQQQrRidqf1IFWK3uNA+w1WshQohBCirQkveg+UWi+l6L3ZtONgZU2XuuwyYyWEEKKNqVz0rrW1TOgvlZ2FzaDdBqvy4nVZChRCCNGWKQXOBOshIavJtdtgJUuBQggh2p3KIcs0rZAVKLOWC0VMtONgJbsChRBCtGM2W6V2OnJGVqy032AluwKFEEIIS/gZWYEyCVkN1H6DVejkdSleF0IIISrY7BKyGqHdBquK4nWpsRJCCCEiCg9Z/uDxDRKyatRug5XUWAkhhBD1UPm0dwlZNWrHwSq0K1CClRBCCFEvErJqJMFKjlsQQgghGi48ZLXzmqz2G6z81lKg065aeCRCCCFEGyGF7+04WBkmbocNpSRYCSGEEDFXJWSZFX0L2/hhpO03WPlNqa8SQgghmkOVc7IqnfjeBtvq1JkslFIOpdQ+pVRe8DGshut+o5RaqZR6JvbDjD1vwJSjFoQQQojmFjrxPTEDkjpBQho43Fa7nTYgmimb4cBrWuvc4GND+AVKqTHAJGAscEwpdXGMxxlz3oAhM1ZCCCFES7LZrL6FiRmQ3NkKWU5Pqw5Z0SSL8cCVSqkVSqnnlVKRlg+nAP/RWmvgU2By+AVKqTuVUvlKqfyCgoLGjToGvAETlwQrIYQQIj6EGkQnpAdDVrr1tmpdP6ujGe1K4GKt9VjACVwR4Zok4GDw9ZNA5/ALtNaztNY5WuucrKysho43ZnwBU45aEEIIIeKRUtbMVUIapHS2ZrRciVZBfJyLpnh9vdY6VMKfDwyIcE0RkBB8PZnoAluL8gakeF0IIYRoFRxu6wEQ8FUUv5tGy44rgmiSxctKqRFKKTvwLWBdhGtWYdVYAYwA9sRmeE3H65caKyGEEKLVcbjAkwrJnSCpY9zVY0WTLB4GXgbWAkuB1Uqp58KuWQSMUko9CdwPvBbTUTYB2RUohBBCtHJ2Z9wFqzqXArXWG7F2Blb2vbBrzOBOwGnAk1rr3bEbYtPwBkxcdpmxEkIIIUTsxOyAUK11KfB2rO7X1HwBA7dTgpUQQgghYqfdJgspXhdCCCFErMVNSxu/38+BAwcoKytrlud7eEo6Hqdi8+bNzfJ8Iv55PB569OiB0+ls6aEIIYRopeImWB04cICUlBR69+7dLI2RzUOnSU900S0toe6LRZuntebEiRMcOHCAPn36tPRwhBBCtFJxsxZWVlZGZmZms4QqsHo+xtlGAtGClFJkZmY224ypEEKItilughXQjKFKY2qNQpKVqNBc339CCCHarrgKVs1FB1/awn6OFhUV8e1vf5tJkybxne98hwcffJC8vLzmHh4Ae/bsifjc9913X/MPRgghhBBRadXBaunSpTz22GMsXbq0Xp9n9YquPkPx9NNPM2DAABYtWoTX6+XNN9+M2Vjrq6Zg9cQTTzT/YIQQQggRlVYbrJYuXcrUqVOZOXMmU6dOrVe4MoNTVuEzVsuXL+eCCy4AYNKkSUybNo3HH3+cKVOmcOONN2IYBlprbrvtNi644AKmTp3K6dOn0Vrz/e9/nylTpjB9+nQMw+pdlJuby1/+8heGD7fOV33llVfKg9Frr73Gn/70JwzD4Oabb2bixIlcffXV+P1+nnzySe677z5efPFFcnNzKSgoKB9jbm5ulTHfc889TJ48mauuuorCwkJefPFFfvKTn3DRRRcxePBgNm7cGHHMQgghhIi9Vhus8vLy8Pl8GIaBz+er15JdcMKqWvH62bNnSUpKAiAxMZEzZ86Qk5PD/Pnz6dChAx988AEnT55k/fr1zJ8/n5kzZ3L69Gnee+89/H4/8+fPJzs7m48++giAw4cPo5Ri/fr1AFx11VV8+eWXAHz66adMnz6dEydOMG3aNObPn09qaiqrV6/m3nvv5YknnmDGjBnk5eWRlZUV8ev48MMPKSsrY+HChVx77bX84Q9/AGDlypV8+umn3H///bz//vsRxyyEEEKI2Gu1wSo3NxeXy4XdbsflclWbyalNTUuBqampFBUVAVBcXExqairjxo0DYPTo0ezcuZPMzExmzJjBZZddxuzZs0lJSWHr1q0sXbqU3NxcFixYwNGjRwHo0KEDP/7xj6vc3+12U1RUxPHjx+nTpw9Op5MPP/yQ6667jl27dlFaWhr117Fp06by8Y0fP778TK6bbroJp9NJdnY2Pp8v4piFEEIIEXutNlhNmDCBefPm8dvf/pZ58+YxYcKEqD/XDL4M/+LHjRtXPvO1cOFC5s6dy6pVqwBYv349vXv3Zv/+/WRmZvLpp5/SvXt35syZw7nnnsuNN95IXl4eTzzxBIMHDwasWS+breqzXHnllTz++ONMnDgRgDlz5jB06FDmzJlD9+7dy69LSEigpKQEqAiC4YYMGcKyZcsAWLZsGUOGDAEon3ULiTRmIYQQQsReqw1WYIWrX/7yl/UKVVDzjNWPfvQjdu3axfnnn09CQgLXXXcdCxcuZMqUKRw9epSrr76aLl268MEHHzBx4kS+/PJLLr74Yr75zW9y6NAhpkyZwoMPPkivXr1qfO6rr76axx9/nOnTpwMwceJE3njjDSZNmsTJkyc5ePAgAKNGjWLr1q1MnjyZN954I+K9pk2bRkJCApMmTeI///kPv/jFLyJeF2nMQgghhIg9VdNsSFPKycnR+fn5Vd63efNmBg0a1CzPX+wNsLOgiL4dk0j2SPsSUaE5vw+FEEK0DkqpVVrrnGiubdUzVg1l1jBjJYQQQgjRGO0yWNW0K1AIIYQQojHaabCykpVNkpUQQgghYqhdBqvQrkCJVUIIIYSIpaiClVKqs1JqTQ0fcyil9iml8oKPYbEdYuzVtCtQCCGEEKIxop2x+j8goYaPDQde01rnBh8bYjO0phOppU14A+ZAINAkz7127VrWrl3bJPcGePHFF3nxxReb9D4vvvgiffr0YcqUKUyaNImvv/46Js8RTYNpaU4thBAintUZrJRSFwHFwJEaLhkPXKmUWqGUel4p5YjlAJtCRfF6RbJqrgbMTR2smssdd9zB/Pnz+f3vf88jjzwSk3tG02BamlMLIYSIZ7WGIKWUC5gJfBt4t4bLVgIXa60PK6X+BVwBvN+YQf3mg6/ZdOhMY25RzeBuqfz6Kutk8oqlwIqPL1++nO9973uA1YB55cqVbNu2Db/fz8KFCzlz5gyffPIJqamp3H777Rw7doxhw4bxzDPPRHy+0tJSrrvuOs6cOUNmZiZvvfUWM2fO5J133gHg5ZdfZt68eZSUlFS730MPPcTy5cspKSkhKyuL119/HYej+l+VYRjcdttt7N27l44dO/L2228DsG7dOi666CKOHDnCm2++yZAhQ7jzzjvZtm0bWVlZ5QeORvrc66+/nlOnTuF0Orn55pvr/HMtLCwkIcGazMzNzeXqq69m9uzZrF+/Hq11tec1TbPG58jNzS0PTVprfvSjH7F27VqcTievv/46b7zxBrNnz+bUqVPk5eXx1ltvlfdRrPy5YDWnXrt2LWlpafzrX//ivffeY926daxbt678z2Xo0KF1fn1CCCFEfdQ1Y3U/8Det9alarlmvtT4cfD0fGBDpIqXUnUqpfKVUfkFBQQOGGjumtgrXK1dYRWrADLBjxw4WLFjANddcw5dffsmsWbMYOmqSmdQAAArLSURBVHQoCxYs4PDhw6xfv54f/OAH5Obmlj8efvhhNm3ahM1mY8GCBXz3u9+lqKiIxx57jPvvv5/777+fefPmAUS8H8DkyZOZP38+nTt35r333uPXv/51lee46667IjZwhupNmCM1iY70uXPmzKFXr1589dVXtZ4eD/D8888zZcoUnn/+eX7zm98A1ZtOR3reaJ/jgw8+IBAIsHjxYn7+85+zatWqmDenFkIIIWKtrmW7i4GLlFJ3AyOVUs9prb8Xds3LSqlHgY3At4DfRbqR1noWMAusk9dre9LQzFJT0WiUUlWWAiM1YAa4/fbbAcobGm/dupUlS5aQl5fHqVOnOHjwIM8++2z159CaoUOHcumllzJgwAAuu+yyiGOJdD+AMWPGADB8+HD27NlTHl4qKyws5MMPP+Ttt9/m2LFj5Q2cKzdh3rVrV5Um0UVFRQwaNKi8+XPlz929ezcjRowAICen9gNm77jjDh588MEq7wtvOh3peU+cOBHVc2zZsoWxY8cCVn9F0zRrvDZceHPqd955h4EDB1b7cxFCCCFirdYZK631BaGidGAt8LhSKryg5mHg5eDHl2qtv2iSkcaQqasfDhregDn0Qz28ofG5557LfffdR15eHo888gjZ2dkRn2PdunVMnDiRzz77jMLCQhYuXAhUb65c0/1WrFgBwJo1a+jfv3/E56ipgXOkMYc3iY70udnZ2eWF6GvWRNwEWqvwptORnjfa5xg4cCArV64E4JVXXmHmzJlAbJtTCyGEELEW9TlWwYC1SWv9YNj7N2qth2uth2mtH4j9EGNPa13tcNBIDZgj+f73v8/HH3/MBRdcwD/+8Q969uwZ8brevXvz1FNPcf7553PkyJHy2ZlLLrmEOXPmMHHiRBYuXFjj/VauXElubi6nTp3iyiuvjPgcNTVwDhepSXSkz7322mvZtm0bubm5bNu2Lao/y9pEet5on+Oqq65CKcUFF1zAyy+/XL7rL5bNqYUQQohYa5dNmPefLKHYF2Bgl9Rmeb76euihh8prqUTzkibMQgghwtWnCXPcH43QFEytUXF87vpDDz3U0kMQQgghRAO0y5Y2Wlc9HFQIIYQQIhbaZbAytZZ2NkIIIYSIuXYZrDTVdwUKIYQQQjRW+wxWmmq7AoUQQgghGqtdBiureL26hx56iEGDBpXvyPvrX//a7GMLaSs9BYUQQoj2pF3uCqxtxuqBBx7g1ltvbeYRVRcKVSNHjmzhkQghhBAiWvEZrO67D2I9WzNyJDzxBGAdEBrNSqDX62XGjBkcOnSIHj16MHv2bFwuF7m5uZx33nmsX7+eTz/9NGIj5bKyMmbMmMGBAwdIS0vjzTffxDRNpk+fTnFxMf3792f27NlRN2sWQgghRPxrp0uBNRevP/roo+VNjv/5z38ydOhQ5s+fz4ABA3jhhRcAq03KhAkT+PTTT4HIjZRnzZrFiBEjWLRoEddeey0bN27k8OHD3HPPPXzxxRfs2bOHo0ePRt2sWQghhBDxLz5nrIIzS01FU72lTUjlpcC77rqLa665BrCa+X788ccADB06tPz9ELmR8pYtW7j22msBmDFjBgB79+7lueeeY/bs2Zw8eZLS0lJGjx4dVbNmIYQQQsQ/mbGqRU3NfJOTk6tcF6mRcuUmwr/73e947rnneP7555k+fTqvvfZaeUPg/9/e3YVIVcZxHP/+0F1HKs0gjdheVgiiiwohSCqZpE26il6gqAQzNej1oosKNQkCCaKboAspoqibuigwjCxDFKkgtwgNS4Je3ChkwqQuRN1/F+dsbdNZPTPOy5lzfh9YmPPMwfOsPx985jxnnn/eYs1mZmZWfJWbWEVEZhHmLGvWrGH//v0sW7aMgwcP/nPnqVlWIeW1a9cyPj5OvV5nfHyclStXMjY2xubNm1m+fDkAExMTuYs1m5mZWfFVrgjz5GSw75c/uGBejYXzal2/ng0WF2E2M7NmrRRhrt4dK5KJpEvamJmZWadVbmI1md6gcxFmMzMz67RCTax6sSw5dQ3fsbJm/pKAmZmdqcJMrGq1Go1Go+v/ufmOlWWJCBqNBrWan7szM7P2FWYfq5GREQ4dOsThw4e7ep3jJyf57egxTjSGmTs8q6vXssFSq9UYGRnpdzfMzGyAFWZiNTQ0xOjoaNev89XPR1j75h5eW3UNSy5f2PXrmZmZWXXkXgqUtEjSl6d4/1VJn0ra0Jmudcex4ycBmDO7MKugZmZmVhKtzC5eAOZmvSHpdmBWRCwFFku6rBOd64ZjJyYBmDPkiZWZmZl1Vq7ZhaTlwF/ArzOcUgfeTl9vB64/4551ydTEaniWn68yMzOzzjrtM1aShoGNwG3AezOcdhYwkb7+HViS8eesA9alh39K+rbl3rbuYuCnrDeufL4HV++N+cAf/e5ED8yYZclUIU9nWS5VyNNZlks7eV6S98Q8D68/BbwcEUdOsffTn/y7THg2GXfCImILsCVvxzpB0uG8W9APKklbImLd6c8cbFXIEqqRp7Mslyrk6SzLpdt55lkKvAl4WNJO4GpJr2Scs5d/l/+uAn7oSO/O3JF+d6AHtva7Az1ShSyhGnk6y3KpQp7Osly6mmdLRZjTydVDwD0RsWFa+zxgN7ADuAW4NiL6fttU0hdVmH1XgbMsD2dZLs6zPJxlZ7T01biIqEfEN9MnVWn7UZIH2D8DbizCpCrV06VH6ypnWR7OslycZ3k4yw5o6Y6VmZmZmc3MmzmZmZmZdYgnVmZmZmYdMtATq7TMzu709RJJH0vaI+mJtO1ZSTvTnwOSnpY0JGlret7q/v4GNiVHlosl7ZD0laSNaZuzLKAcWWa1OcuCkTRf0geStkt6V9JwVumyvG3WPy1k+c/YTY89LtswsBMrSQuA10k2JwV4CbifZNuHOySNRsSm9IH7OrAPeAN4FNgbEdcBd0o6p/e9t+nyZAk8AjwTEVcDKySdj7MsnJxZZrU5y+K5F3gxIm4mqbpxN02ly7LKmQ1SibMKyZNl89gFj8u2DOzECjgJ3AUcTY/Pi4ifI3kavwHMmzpR0jXAoYiY4L/ld3YB/mpp/+XJsgFcKWkRMIdkv5U6zrJo8mSZ1VbHWRZKRLwcER+lh+cD9/H/0mX1nG3WRzmzbB674HHZljw7rxdSusUD03aD3yPpEZKSOpcCX087/XFgU/q6ufzOom731U4tZ5azgceAEeAT4ATOsnByZpnV5iwLStJSYAHJxs/NpcuyypmdtsSZ9cepsswYu+Bx2ZZBvmPV7EHgAMmS0fPpp2EknQssjIjv0/NOW37H+i4ry6eAVRGxniS/MZzlIMjKMqvNWRaQpPNIlm5Xk51R3jbrsxxZZnGWbSjNX1JEnASmCju/Ne2tW4Ft046LWn7HUjNkOQpcJKlG8gk4cJaFl5XlDPk6y4KRNAy8AzwdET+SnVHeNuujnFlmcZZtGNilwBk8Bzw5dbcqtQJ4Ydrx68A2STcAVwCf97B/ll9zlpuAnSTPB7xPshz4Hc5yEGSNy+Y2j8vieYDkQ8x6SeuB14CVki4kLV1G8gFnd4426688WWbxuGxDJXdeT/8xXQ98WKDyO9YGZ1kezrL40m+OjQG7IuLXVtqsWPJm5HHZukpOrMzMzMy6oTTPWJmZmZn1mydWZmZmZh3iiZWZmZlZh3hiZWZmZtYhnliZmZmZdcjfkZUtqI/z92AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,4))\n",
    "\n",
    "# Plot the results\n",
    "df['lff'].plot(ax=ax, style='k.', label='Observations')\n",
    "predict.predicted_mean.plot(ax=ax, label='One-step-ahead Prediction')\n",
    "predict_ci = predict.conf_int(alpha=0.05)\n",
    "predict_index = np.arange(len(predict_ci))\n",
    "ax.fill_between(predict_index[2:], predict_ci.iloc[2:, 0], predict_ci.iloc[2:, 1], alpha=0.1)\n",
    "\n",
    "forecast.predicted_mean.plot(ax=ax, style='r', label='Forecast')\n",
    "forecast_ci = forecast.conf_int()\n",
    "forecast_index = np.arange(len(predict_ci), len(predict_ci) + len(forecast_ci))\n",
    "ax.fill_between(forecast_index, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], alpha=0.1)\n",
    "\n",
    "# Cleanup the image\n",
    "ax.set_ylim((4, 8));\n",
    "legend = ax.legend(loc='lower left');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "    Commandeur, Jacques J. F., and Siem Jan Koopman. 2007.\n",
    "    An Introduction to State Space Time Series Analysis.\n",
    "    Oxford ; New York: Oxford University Press.\n",
    "\n",
    "    Durbin, James, and Siem Jan Koopman. 2012.\n",
    "    Time Series Analysis by State Space Methods: Second Edition.\n",
    "    Oxford University Press."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
